March  31, 1993 


Dr.  Edwin  P.  Rood 
Scientific  Officer  Code:  1132F 
Office  of  Naval  Research 
800  North  Quincy  Street 
Arlington,  VA  22217-5000 

Dear  Dr.  Rood: 

The  purpose  of  this  letter  is  to  transmit  the  ninth  quarterly  report  for  ONR  Grant 
N00014-91-J-1271,  "An  Experimental  Study  of  a  Plunging  Liquid  Jet  Induced  Air 
Carryunder  and  Dispersion"  (Lahey  &  Drew  -  Co-PI). 

During  this  report  period  a  technical  paper  entitled,  "A  Numerical  Simulation  of 
Two-Phase  Jet  Spreading  Using  an  Eulerian-Lagrangian  Technique,"  was 
accepted  for  presentation  and  publication  at  the  1993  Winter  Annual  Meeting  of 
the  ASME.  A  draft  of  this  paper  was  previously  transmitted  to  you  by  the  eighth 
quarterly  report  and  it  has  only  needed  to  be  slightly  modified  to  resolve  reviewer 
comments. 


This  report  period  was  primarily  concerned  with  performing  mechanistic  CFD 
predictions  of  a  spreading  two-phase  jet  using  a  multidimensional  two-fluid 
model,  and  the  comparison  of  these  predictions  with  the  data  we  have  previously 
acquired.  The  attached  paper,  entitled,  "A  Numerical  Simulation  of  a  Turbulent 
Two-Phase  Jet  Using  a  Multidimensional  Two-Fluid  Model"  has  been  prepared 
for  submission  to  the  Int  </.  Numerical  Methods  in  Fluids.  As  you  can  see  the 
two-fluid  model  predictions  agree  fairly  well  with  the  data,  however  some  trends 
have  not  yet  been  fully  predicted.  Moreover,  more  work  is  still  required  to  properly 
model  the  air  capture  and  spreading  phenomenon  as  the  liquid  jet  enters  the  pool. 
Nevertheless,  these  results  are  quite  exciting  since  they  show  that  mechanistic 
predictions  of  two-phase  flow  phenomena  in  free  jets  are  indeed  possible. 

Finally,  as  I  have  previously  reported  to  you,  in  early  March  I  met  with  Professors 
Lasheras  (UC-SD)  and  Loth  (U.  Illinois)  in  Pasadena  to  discuss  research 
collaboration.  In  addition,  I  visited  Professor  Lasheras'  laboratory  in  San  Diego. 
It  appears  that  research  collaboration  would  be  mutually  beneficial  and  Professor 
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Drew  and  I  plan  to  formally  establish  synergistic  research  collaboration  in  a 
future  extension  of  this  research  project. 


If  you  need  any  further  information  concerning  this  report  please  don’t  hesitate  to 
contact  me  [(518)  276-8579]  or  Professor  Drew  [(518)  276-6903]. 

Sincerely  yours. 


Dr.  R.T.  Lahey,  Jr.  / 

The  Edward  E.  Hood,  Jr.  Professor  of  Engineering 
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INTRODUCTION 

A  good  understanding  of  the  air  carryxuider  and  bubble  dispersion  process 
associated  with  a  plunging  liquid  jet  is  vital  if  one  is  to  be  able  to  quantify  such 
diverse  phenomena  as  sea  surface  chemistry,  the  meteorological  significance  of 
breaking  ocean  waves  (eg,  mitigation  of  the  "greenhouse"  effect  due  to  the 
absorption  of  CO2  by  the  oceans),  the  performance  of  certain  type  of  chemical 

reactors,  and  a  number  of  other  important  maritime-related  applications. 

The  absorption  of  greenhouse  gases  into  the  ocean  has  been  hypothesized  to 
be  highly  dependent  upon  the  air  carryunder  that  occurs  due  to  breaking  waves. 
This  process  can  be  approximated  with  a  plunging  liquid  jet  [Monahan,  1991; 
Kerman,  1984].  Moreover,  the  air  entrainment  process  due  to  the  breaking  bow 
waves  of  surface  ships  may  cause  long  (ie,  up  to  5  km  in  length)  wakes. 
Naturally,  easily  detectable  wakes  are  undesirable  for  naval  warships.  In 
addition,  the  air  carryrmder  that  occurs  at  most  hydraulic  structures  in  rivers  is 
primarily  responsible  for  the  large  air/water  mass  transfer  that  is  associated  with 
these  structures  [Avery  and  Novak,  1978],  Also,  air  entrainment  plays  an 
important  role  in  the  slug  flow  regime.  In  particular,  the  liquid  film 
surrounding  a  Taylor  bubble  has  a  flow  in  the  opposite  direction  from  the  Taylor 
bubble.  This  liquid  film  can  be  thought  of  as  a  plimging  liquid  jet  that  produces  a 
surface  depression  in  the  rear  part  of  the  Taylor  bubble.  When  the  annular  liquid 
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jet  exceeds  a  critical  velocity,  it  entrains  small  bubbles  from  the  air  in  the  Taylor 
bubble.  These  entrained  bubbles  follow  the  Taylor  bubbles  in  the  liquid  slug. 

Single-phase  turbulent  jets  represent  an  important  class  of  free  shear  flows 
that  have  been  studied  in  the  past  to  develop  and  test  turbulence  models 
[Abramovich,  1963].  More  recently,  turbulent  jets  have  been  evaluated 
numerically  using  computational  fluid  d3mamic  (CFD)  techniques  for  various 
turbulence  models.  Rodi  [1984]  presented  results  using  the  classical  k-e  model  of 
Gibson  &  Launder  [1976],  and  showed  that  k-e  models  may  not  accurately  predict 
jet  spreading.  Rodi  [1984]  proposed  that  the  constant  in  the  model  for  turbtdent 

viscosity  was  really  a  function  of  the  ratio  between  turbulent  production  and 
dissipation.  This  involved  the  development  of  a  new  function  which  produced 
better  results.  Sini  and  Dekeyser  [1987]  solved  the  single-phase  turbulent  jet  using 
Rodi's  k-e  model  [Rodi,  1984].  This  model  compared  favorably  with  the 
experimental  results  of  the  single-phase  turbulent  jet  as  well  as  with  other  more 
detailed  algebraic  stress  models.  Hence  it  appears  that  in  some  cases  turbulent 
nonisotropy  is  not  important  and  need  not  be  modeled. 

Significantly,  it  has  been  found  that  single-phase  turbulent  jet  data  can  be 
used  for  the  assessment  of  turbulence  models  because  one  does  not  have  to 
constitute  complicated  turbulent  closure  laws  near  solid  (no  slip)  boundaries. 
Indeed,  due  to  the  absence  of  walls  and  the  associated  shear  boundary  conditions 
the  turbulent  jet  is  probably  the  simplest  non-trivial  case  to  analyze. 
Interestingly,  the  same  conclusions  can  be  reached  for  a  two-phase  turbulent  jet. 

In  most  of  the  previously  mentioned  research,  the  flow  field  was  considered 
to  be  thin  in  the  lateral  direction  and  the  flows  are  characterized  by  a  relatively 
small  lateral  velocity  when  compared  to  the  streamwise  velocity.  This  is 
equivalent  to  considering  the  flow  (ie,  the  liquid  jet)  to  be  a  boundary  layer.  Hence, 
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the  flow  is  often  assumed  to  have  a  negligible  pressure  gradient  in  the  lateral 
direction,  and  the  pressure  in  the  boundary  layer  can  be  imposed  by  the  external 
flow.  From  the  mathematical  point  of  view,  this  approximation  changes  the 
problem  from  an  elliptic  to  parabolic  one  in  the  streamwise  direction.  From  the 
computational  point  of  view,  when  using  a  boimdary  layer  approximation  the 
pressure  distribution  has  to  be  prescribed.  That  is,  the  static  pressure  as  a 
function  of  the  streamwise  coordinate  has  to  be  known  and  supplied  to  the 
computational  fluid  dynamic  (CFD)  code. 

Depending  on  the  variables  of  interest,  approximating  the  jet  as  a  boundary 
layer  may  be  a  reasonable  approximation.  However,  the  lateral  velocity  at  the 
edge  of  the  jet  (i.e.,  the  entrainment  velocity)  may  be  much  larger  than  the 
streamwise  velocity.  Moreover,  for  a  planar  jet  the  entrainment  velocity  remains 
finite  as  the  integration  domain  in  the  lateral  direction  is  enlarged. 

In  this  work  we  did  not  make  the  boundary  layer  approximation.  Rather, 
we  assumed  an  elliptic  problem  for  both  the  gas  and  the  liquid  phases.  Solving 
the  partial  differential  equations  as  an  elliptic  system  increases  the  complexity  of 
the  problem  but  provides  more  detailed  and  accurate  information  on  the  flow  field 
than  a  parabolic  scheme.  One  of  the  advantages  of  the  elliptic  solution  is  that  the 
streamwise  pressure  distribution  does  not  have  to  be  provided.  This  has  two 
effects,  first  we  are  able  to  compute  recirculation,  and  second,  we  are  able  to 
calculate  the  buoyancy-induced  reversal  of  the  bubbles.  These  two  important 
effects  cannot  be  computed  using  a  parabolic  approach. 

In  parabolic  single-phase  jet  calculations  using  a  k-e  model,  Sini  & 
Dekeyser  [1987],  and  Hossain  &  Rodi  [1982]  found  satisfactory  agreement  between 
calculations  and  experiments,  except  for  a  small  region  near  jet  inlet.  Their  good 
results  are  due  in  part  to  the  fact  that  they  analyzed  a  free  jet  which  was  only 
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weakly  nonisotropic.  The  next  level  of  complexity  for  simulating  turbulent  two- 
phase  flows  is  the  use  of  an  Algebraic  Stress  Model  (ASM).  Several  different 
models  have  been  proposed  in  the  literature.  For  the  particular  case  of  a  planar 
jet,  performance  of  the  ASM  is  similar  to  the  k-e  model.  One  might  expect  that 
the  nonisotropic  ASM  model  would  produce  better  results  than  the  isotropic  k-e 
model  for  the  pressure  distribution.  However,  when  the  ASM  proposed  by  Gibson 
&  Launder  [1976]  was  used  for  the  evaluation  of  a  planar  jet  the  computed 
streamwise-pressure  distribution  [Bergstrom,  1992]  was  virtually  the  same  as 
when  a  k-e  model  was  used. 

Thus,  for  simplicity,  we  have  used  a  k-e  model  in  the  computations 
presented  in  this  paper.  We  note  that  the  turbulence  present  in  the  liquid  has  two 
components  in  a  two-phase  jet.  One  component,  the  shear-induced  turbulence,  is 
due  to  viscosity  and  it  is  present  in  both  single  and  two-phase  flows.  The  other 
component  is  the  bubble-induced  turbulence  due  to  slip  between  the  bubbles  and 
the  surrounding  liquid,  and  it  only  occurs  in  two-phase  flows. 

THE  TWO-FLUID  MODEL 

Different  phenomenologically-based  models  for  two-phase  flows  have  been 
proposed  in  the  past.  The  drawback  of  many  of  these  models  is  that  they  are  only 
applicable  to  particular  problems.  On  the  other  hand,  mechanistically-based 
models,  commonly  known  as  Two-Fluid  Models  (TFM),  have  been  proposed  [Ishii, 
1975;  Delhaye,  1968  and  Drew  &  Lahey,  1979]. 

For  an  adiabatic  plunging  liquid  jet  entraining  air  bubbles,  we  have  the 
following  local,  instantaneous  conservation  equations: 
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Mass. 

dpi, 

-^  +  V*(pj,xi,)  =  0  (k  =  g,l)  (1) 

Momentum 

^(pi,Xk)  +  V*(pkYkEk)  =  V*Sc  +  Pk£  {k  =  g,/)  (2) 

where, 

Sc^-Pkl+Jk  ^3) 

and,  pg,  p^,  Xg,  S!/.  Pg.  P/»  Jg.  phasic  densities,  velocities,  pressures  and 

shear  stresses  of  the  gas  and  liquid  phases,  respectively. 

These  local,  instantaneous  conservation  equations  may  be  appropriately 
averaged  to  obtain  the  two-fluid  model.  Ishii  [1975]  and  Delhaye  [1968]  have 
proposed  time  and  spatial  averages.  However,  the  ensemble  average  [Drew  & 
Lahey,  1979]  seems  to  be  the  most  fundamental  type  of  averaging.  An  ensemble 
average  is  defined  as, 

f(i;t)  =  )f(x,t;|)P(^)d5  (4) 

3 

were  f  is  the  function  we  want  to  average,  x  is  the  spatial  coordinate,  t  is  time,  %  is 
a  parameter  that  determines  a  particular  realization,  P{4)  is  the  probability 
density  function,  and  S  is  the  set  of  all  realizations.  We  note  that  P(^)  satisfies, 

|f(x,t;?)d5=1.0  (5) 

S 

Let  us  define  the  phase  indicator  function,  X]((x,t),  such  that  it  is  unity  if  phase-k 
is  present  at  2  and  time  t,  and  is  zero  otherwise. 
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Thus,  the  volume  fraction  of  phase-k,  a^,  is, 

Ok  =  (6) 

The  physical  interpretation  of  aj^  is  as  follows;  ajj(2{jt)  is  the  fraction  of  all 

realizations  in  which  phase-k  is  present  at  location  x  at  time  t. 

We  refer  the  reader  to  Arnold  [1988]  and  Park  [1992]  for  complete  details  on 
the  derivation  of  the  two-flmd  model.  We  present  here  only  the  final  restilts. 

The  ensemble-averaged  continuity  equation  for  phase-k  is: 


3(akPk) 


at 


+  ^*(okPkVk)  =  ') 


(k  =  g  or  t)  (7) 


where, 


Pk 


^kPk 


(8) 


is  the  ensemble-averaged  phasic  density. 


In  this  work  we  assume  that  the  gas  and  liquid  densities  are  both  constant,  thus. 


3^ 

Pk  =  Pkak  =  Pk 


(9) 


Next,  Yk*  ensemble-averaged  phasic  velocity,  is  given  by: 

“kPk  “k 

The  ensemble-averaged  momentum  equation  for  phase-k  is: 
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+“kPki'^— ki 


(k  =  g  or  /) 


where,  Xj,.,  g,  Mki  ®re  the  averaged  stress  tensor,  the  acceleration  of  gravity  and 
the  total  interfacial  forces,  respectively.  They  are  given  by: 

f  (12) 

«k 

T^  =  ”PkSiik 

where,  is  the  fluctuation  in  the  velocity  of  phase-k,  and, 

The  ensemble  averaged  interfacial  jump  conditions  are: 


We  note  that,  for  monodispersed  spherical  bubbles,  Laplace's  equation  yields, 

i~-ij:  =  2a/Rb  (16) 

where,  a  is  the  surface  tension  and  the  mean  bubble  radius.  The  bubble's 
surface  stress  tensor,  is  given  by  [Park,  1992], 
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where,  =  Vg  -  v^,  is  the  average  relative  velocity,  and,  using  inviscid  flow 
theory,  the  coefficients  a^  and  in  Eq.  (17)  can  be  analytically  computed  for 

spherical  bubbles  to  be: 


^,4 


(18) 


TURBULENCE  MODELING 

The  total  Reynolds  stress  tensor  for  the  continuous  liquid  phase  is  given  by 
superposition  as. 


Re  Re  Re 
It  =J/(BI)  +I#(SI) 


where,  for  bubbly  two-phase  flows  the  bubble-induced  shear  stress  is  given  by, 
Re 


(19) 


(20) 


The  bubble-induced  turbulent  kinetic  energy,  k^(Bi),is  given  by  Eq.  (23a),  where  for 
potential  flow  around  a  sphere, 

f  9/10  0  0  ^ 

(21) 


A(bi)  = 


0  9/10  0 

^  0  0  12/10  J 

Re 


We  note  that  J^(si)  is  the  shear-induced  Re3niolds  stress  which  comes  from  the 
classical  k-e  model  [Rodi,  1984]: 

^  V.[a,«Jvk«s„]  +  a/P,  -  e,) 


Dt 


r 

r  T 1 

=  V* 

«/ 

PrJ 

Ve, 

L  t  j 

> 

+  a/P^-ep 


(22a) 


(22b) 
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where, 

1  —2 
^/(BI)  =  2“g 

(23a) 

=  ^£(SI)  ^ 

(23b) 

2 

p  p  MSI) 

(23c) 

=  "0^  (VY/  +Y/  V):Vy^ 

(23d) 

e/P/ 

^  Msi) 

(23e) 

h  -  ^2e 

(23f) 

PrJ  =  1.3 

(23g) 

and  [Rodi,  1984],  =  0.5478,  =  0.1643,  =  1.44,  €2^  = 

1.92  are  the  single- 

phase  flow  k-e  model  coefficients. 

Using  these  results,  the  Reynolds  stress  tensor  is  given  by; 


After  averaging,  we  have  more  unknowns  than  equations.  Thus,  we  need 
to  constitute  some  of  the  averaged  quantities  in  terms  of  the  state  variables, 
and  pjj.  This  closure  process  is  necessary  because  we  have  lost  information  due 

to  averaging  and  we  must  reintroduce  the  essential  physics  which  was  lost. 

Cell  model  averaging  was  successfully  applied  to  a  dilute  mixture  of  liquid 
and  gas  bubbles  by  Arnold  [1988]  and  Park  [1992].  As  can  be  seen  in  Figure-1,  we 
divide  the  flow  field  into  "cells",  each  of  which  have  only  one  spherical  bubble 
inside.  Using  inviscid  flow  theory  we  may  compute  the  pressure  distribution 
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around  the  bubbles  and  thus  deduce  the  various  interfacial  forces.  The 

assumptions  are  that  both  phases  are  invisdd,  incompressible  and  have  constant 

thermophysical  properties,  the  bubbles  are  spherical  and  can  be  treated  as  a 

dilute  dispersion  of  spheres,  the  non-uniformities  in  the  distribution  of  the 

dispersed  phase  are  small  and  the  velocity  gradients  of  both  phases  are  small. 

Figure- 1  shows  a  typical  cell  for  one  particular  realization  of  the  ensemble. 

Note  that  x  is  the  vectoral  location  under  consideration,  z.  is  the  position  of  the 

center  of  the  bubble  in  a  particiilar  realization  of  the  ensemble,  and  x'  =  x  -  i  is  the 

position  of  the  bubble  with  respect  to  x,  the  location  under  consideration.  We  need 

-Re  -  - 

to  constitute  the  following  quantities:  2k  >  St'  M.ki»  Pki  terms  of  the  state 
variables. 


Using  the  results  of  Arnold  [1988]  and  Park  [1992]  the  resultant  two-fluid  phasic 
momentum  equations  are: 

Momentum  Conservation  -  Gas  Phase 


at  '“gPe5:g)  +  ^  *  <“gPg!^g%>  =  ■  “g^' 


^vm®gPt 


-  Crot^gP/^r  X  V  X  Vg  -  X  V  x 


-  (Cj  +  C2  2Cp  -  2b5)  ttgp^r  •  J  -  Cg)  Ogp^r  • 

—  _  Cjj  _  — 

+  (ag  -  C2)agp/V  •  y^)  y^  +  ttgPgg  - p^  Aj'^y,  \Xr^  -  C-rDP/k^^ag  (25a) 


u 


Momentum  Conservation  -  Liquid  Phase 

a  -  -  -  -  -  2 

5t  (“/P/Xp  +  V  •  +  (Cp  +  bg  +  bp  pP  xp 

+  C™«gP'[(s  +  ^g*'')^g-(s  +  ^< 

+  C^ot-agP^r  X  ^  X  ig  +  CLP^gXr  X  V  X  X/  +  (Cg  +  atkLgPf  Xr  (V  •  Xf) 

+  (Cl  +  Cg  +  2bp  OgP^Xr  •  VxJ  +  (Cg  +  ap  Cgp^  Xy*  VXf 

—  —  —  — 

+  (ag  +  ap  p/  ( Xr  •  Vttg)  Xr  +  a^P/  g  +  X  '  +  CxDP/k^Vag  (25b) 


where  both  phases  were  assumed  to  be  incompressible. 

For  numerical  purposes  it  is  important  to  rewrite  Eq.  (25a)  as  follows.  For 
the  air/water  flows  imder  consideration  Eq.  (7)  can  be  employed  to  show  that  the 
left  hand  side  of  the  gas  phase  momentum  equation  can  be  rewritten  in 
Lagrangian  form  as: 


^  («gPg^g)  + 


Dv 


^“gPg'^g^g^  ~  Pg®g  Dt 


K 


(26) 


Grouping  the  right  hand  side  of  Eq.  (26)  with  one  part  of  the  virtual  mass  force  in 
Eq.  (25a)  we  obtain: 


Dv  Dv  Dv 

Pg“g "Dt  +  “g  m  =  ‘Pg ^vmP/) “goT 


(27) 


Notice  that  we  have  kept  the  other  part  of  the  virtual  mass  force  associated  with 
liquid  phase  acceleration  in  tlie  second  term  on  the  right  hand  side  of  Eq.  (25a).  It 
has  been  found  that  writing  the  gas  phase’s  acceleration  as  in  Eq.  (27)  promotes 
numerical  convergence. 
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Equations  (25)  were  numerically  evaluated  using  the  finite  difference 
formulation  of  Patankar  [1980]  in  the  well-known  PHOENICS  code.  First  the 
domain  of  interest  (DOI)  was  subdivided  into  a  Cartesian  grid,  as  shown  in  Pig.  2. 
The  dependent  variables  were  calculated  and  stored  at  discrete  points  on  the  grid. 
To  prevent  pressure  "checker  boarding"  [Patankar,  1980]  a  staggered  grid  was 
used.  The  velocities  are  calculated  at  the  locations,  shown  by  the  arrows  in  Fig.  2 
that  are  between  pressure  points  (P).  The  cell  surroxmding  point  P  is  often  called 
the  continuity  cell.  The  velocities,  Vy  and  v^,  were  computed  at  the  arrow 

locations  and  the  pressure  and  void  fraction,  p^  and  Ug,  were  computed  at  the 
continuity  cell  center  (P). 

According  to  the  differencing  procedure,  the  conservation  equations  were 
first  integrated  over  the  control  volume  that  surrounds  the  node.  The  resulting 
integrals  were  then  approximated  using  the  nodal  values  and  algebraic  difference 
equations  were  obtained,  where  the  discrete  equations  had  an  implicit 
formulation. 

Implementation  of  the  various  source  terms  required  special  attention.  We 
have  used  lagged  quantities  for  them,  (i.e.,  values  from  the  previous  iteration). 
Even  though  there  are  many  terms  on  the  right  hand  side  of  the  conservation 
equations,  all  of  them  can  be  written  as  a  function  of  the  following  quantities:  ttg, 

Vttg,  2}^  and  Vyjj.  Notice  that  is  naturally  given  at  the  continuity  cell  center 
and  Vttg  is  computed  using  a  first  order  Taylor  series  expansion  at  locations  w,  e, 

s  and  n  (Fig.  2).  The  velocity  component,  vj^y,  was  given  at  locations  w,e  and 
was  given  at  locations  n,s.  Finally,  Vv|^  =  Vj^j  j  was  computed  using  a  first  order 

Taylor  series  expansion  and  the  locations  are  P,  nw,  ne,  sw,  se  depending  on  the 
velocity  components  j  and  the  coordinate  i.  Because  these  four  basic  quantities  are 
naturally  given  at  different  locations,  we  need  to  perform  some  kind  of  averaging 
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before  computing  the  source  terms.  To  illustrate  this  let  us  take  as  an  example 
the  lift  force  in  the  y-direction.  The  lateral  force  is  proportional  to: 


were  Og  is  given  at  the  continuity  cell  center,  P  (Fig.  3),  VOg  at  w,  e,  s  and  n,  is 
given  a  w,e,  v^  ^  given  at  nw  and  v^^.y  *®  given  at  ne.  It  is  not  consistent  to 

multiply  quantities  given  at  different  locations.  After  some  experimentation  we 
obtain  acceptable  convergence  speed  and  accuracy  with  the  following  criteria.  All 
tensor  components  with  repeated  indices  (Tjj,  T22.  •••)  were  computed  at  the 

continuity  cell  center  P.  Off-diagonal  tensor  components  (T12.  •••)  were  computed 


at  corresponding  continuity  cell  corners  (eg,  ne  for  the  case  of  index  12). 

In  order  to  numerically  evaluate  the  two-fluid  model  we  also  needed  the 
appropriate  boundary  and  initial  conditions.  There  is  no  general  theory  for  the 
type  of  boundary  and  initial  conditions  that  a  system  of  nonlinear  partial 
differential  equations  requires  in  order  to  have  a  unique  solution.  Many 
researchers  in  the  past  have  evaluated  two-phase  models  as  initial  value 
problems  using  parabolic  numerical  techniques,  thus  it  was  not  necessary  to 
specify  boundary  conditions  at  the  outlet  of  the  integration  domain.  However,  as 
discussed  previously,  when  this  approach  is  used,  one  has  to  specify  the  pressure 
distribution  in  the  integration  domain.  For  many  cases  a  hydrostatic  pressure 
distribution  was  a  good  enough  approximation.  However,  when  a  parabolic 
scheme  is  used,  one  caimot  compute  flow  recirculation  nor  buoyancy-induced  gas 
reversal,  an  important  feature  of  the  two-phase  jets  which  are  of  interest  in  this 
study. 

As  noted  previously,  in  this  work  we  used  an  elliptic  (ie,  boundary  value) 
calculational  scheme.  That  is,  we  numerically  evaluated  the  full  two-fluid  model 
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and  the  associated  k'E  model  using  appropriate  boundary  conditions  at  the  inlet 
and  the  exit  of  the  flow  domain.  This  complicated  the  evaluation  procedure, 
however  only  in  this  way  could  we  obtain  an  accurate  prediction  for  the  gas  phase. 

Figure-3  shows  the  integration  domain.  Note  that  we  have  refined  the  grid 
near  the  symmetry  plane  of  the  planar  jet  because  the  gradients  are  the  steepest 
there.  In  the  axial  (i.e.,  z)  direction  we  have  increasingly  larger  cells  because  of 
the  decreasing  gradients.  The  boundary  conditions  which  have  been  used  are: 

mLEI(z  =  0,-|^y^|> 


Gas  mass  flux  =  a  Pg  V 
Liquid  mass  flux  =  (1  -  a)  V 

5^  =  0 


V  =  v 

3-2 

Kinetic  energy  s 

-  3 

Dissipation  s  e  =  ~yl —  =  1.8256  /h 

C  h 

4 


OUTLET  (z  =  Z) 


p  =  p/g  Z  cos  8 


du 


^X._ 


dz 


=  0 


du 


dz 


du 


Iz 


dz 


=  0 


du 


dz 


(29a) 

(29b) 

(29c) 

(29d) 

(29e) 

(29f) 

(29g) 

(29h) 


(30a) 

(30b) 

(30c) 
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i=0  if=0  (30d) 

It  should  be  noted  that  the  botindary  conditions  associated  with  the  mass  balances 
are  only  required  at  the  inlet.  This  can  be  understood  based  on  the  fact  that  the 
mass  balance  is  a  first  order  partial  dififerential  equation.  Moreover,  given  the 
velocity  field,  the  equation  may  be  solved  for  the  void  fraction  along  the  model's 
characteristics.  With  the  void  fraction  given  at  the  inlet,  the  void  fraction  field  can 
be  readily  evaluated. 

Because  the  velocity  fields,  as  well  as  the  turbulent  kinetic  energy  and  the 
dissipation,  are  not  well  known  at  the  outlet,  we  specified  natural  boundary 
conditions  there.  That  is,  we  set  the  gradients  of  these  variables  in  the  flow 
direction  equal  to  zero.  Fortunately,  in  our  case,  the  outlet  boundary  conditions 
were  found  to  have  a  negligible  effect,  except  for  the  last  two  rows  of  cells. 

The  numerical  residts  presented  herein  correspond  to  a  planar  liquid  jet 
impacting  a  liquid  pool  using  the  two-fluid  model  previously  discussed.  The 
liquid  jet  velocity  at  the  location  of  impact  was  u^^  ^  ^  constant  bubble 

diameter  of  2mm  was  assumed.  The  initial  velocity  profile  was  uniform.  The  jet 
width,  h,  was  4.03  mm,  and  the  void  fraction  at  the  location  of  impact  was 
assumed  to  be  uniform  and  equal  to  5%.  The  inlet  turbulent  intensity  was  3%  and 
the  inlet  dissipation  was  computed.  The  inclination  angle  of  the  jet  was 
measured  with  respect  to  the  vertical  plane  (i.e,  6  =  0“  means  a  vertical  planar 
jet).  The  integration  domain  had  an  extent  of  y  =  0.2  m  in  the  lateral  direction  and 
z  =  0.25  m  in  the  axial  direction.  The  k-e  model  for  turbulence  employed  the 
constant  values  suggested  by  Rodi  [1984]  for  single-phase  flow. 

We  have  presented  results  corresponding  to  a  vertical  liquid  jet  (0  =  0®,  z  = 
0.225  m;  note,  y  =  0.1  m  is  the  jet's  plane  of  symmetry).  Both  velocity  profiles  show 


a  Gaussian-like  profile.  We  see  that,  due  to  buoyancy,  the  liquid  velocity  is  always 
higher  than  the  gas  velocity,  a  well-known  characteristic  of  bubbly  two-phase 
downflows.  Moreover,  the  relative  velocity  was  approximately  0.3  m/s  which  is 
very  dose  to  the  terminal  rise  veloldty  of  the  single  bubbles. 

Figure-5  shows  the  axial  liquid  velocity  as  a  function  of  lateral  position  for 
different  axial  positions.  The  curve  labeled,  z  =  0.00125  m,  is  right  under  the 
location  of  jet  impact  and  one  can  see  an  almost  uniform  velocity  profile  (u^  = 

5  m/s).  As  we  move  down  in  the  pool  the  jet  is  dispersed  due  to  momentum 
interchange  with  the  surroimding  fluid.  The  curve  z  =  0.225  m  is  the  same  one 
shown  in  Figure  4  for  liquid  velocity. 

Figure-6  shows  the  turbulent  kinetic  energy,  k^,  as  a  function  of  lateral 

position  for  z  =  0.225  m.  The  curve  shows  the  characteristic  relative  minimum  in 
k^  at  the  symmetry  plane. 

In  Figure-7  the  liquid  velodty  field  has  been  plotted  as  a  function  of  axial 
and  lateral  position.  The  length  of  the  arrows  is  proportional  to  the  liqmd  velodty 
at  the  location  of  the  arrow's  tail.  The  arrow’s  tail  is  located  at  the  center  of  the 
computational  cell.  The  arrow  scale  is  in  the  lower  left  comer  (u^/  =  2  m/s).  The 

spreading  of  the  jet  can  be  easily  seen  in  this  plot,  as  noted  previously.  Near  the 
location  of  jet  impact  (z  =  0)  the  axial  velodty  is  almost  uniform.  Because  of  the 
momentum  interchange  between  the  jet  and  the  surrounding  fluid,  liquid  is 
entrained  in  the  lateral,  y,  direction.  Finally,  one  may  note  the  formation  of  two 
weak  redrculation  zones  near  the  y-boxmdaries  for  large  z. 

Figure-8  shows  a  contour  plot  of  the  axial  liquid  velodty.  The  lines  connect 
positions  with  the  same  axial  velodty  (equivelodty  lines).  The  outer  curve 
corresponds  to  =  0.25  m/s,  and  the  difference  between  successive  lines  is  0.25 
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Figure-9  shows  a  contour  plot  of  the  void  fraction  for  one  half  of  the  jet.  The 
outer  line  connects  points  with  the  void  fraction  a  =  0.25%. 

Figure- 10  shows  a  vector  plot  of  the  liquid  velocity  field  for  an  inclined 
planar  jet  (0  =  45°).  We  have  rotated  the  integration  domain  45“  in  order  to  have 
the  plane  y  =  0  aligned  with  the  jet  orientation.  This  was  done  to  minimize 
numerical  diffusion.  It  can  be  seen  that,  as  expected,  the  rising  gas  drags  the 
liquid  away  from  the  centerplane. 

Figure-11  shows  a  vector  plot  of  the  gas  velocity  field  for  an  inclined  planar 
jet  (6  =  45“).  Of  particular  interest  are  the  results  shown  in  the  upper  right 
comer,  the  gas  velocity  (weighted  by  the  local  void  fraction)  at  the  y-boundary 
which  shows  flow  reversal  of  the  gas. 


COMPARISON  WITH  EXPl 


NTS 


The  results  presented  in  the  previous  section  have  inlet  boundary 
conditions  that  are  the  simplest  to  implement  for  the  two-phase  jet.  The  average 
axial  liquid  velocity,  and  the  liquid  turbulent  kinetic  energy,  k^,  were 

assumed  to  be  uniform  in  the  liquid  jet  cross  section.  These  are  good 
approximations  for  the  experiments  we  have  compared  the  calculation  against 
[Bonetto  &  Lahey,  1993].  However,  the  assumption  of  a  uniform  inlet  void  fraction 
distribution  was  too  crude.  We  re-evaluated  the  two-fluid  model,  but  instead  of 
using  the  boundary  conditions  given  in  £q.  (29),  we  used  the  average  axial  liquid 
and  gas  velocities,  v^^  liquid  turbulent  kinetic  energy,  and  the 

gas  void  faction  that  were  actually  measured  at  the  inlet  of  the  integration 
domain.  Figure  12  shows  the  average  axial  Hquid  velocity,  v^^*  open  circles 
are  the  experimental  values  at  the  inlet  of  the  domain.  The  solid  curve  is  the 
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computed  v^2  axial  position  z  =  2.5  mm  (i.e.,  at  the  first  velocity  node),  and 

y  =  0.1  m  corresponds  to  the  planar  jet's  centerplane. 

Figure  13  shows  the  computed  mm.  The  open  circles  are 

experimental  points.  We  can  see  that  the  agreement  is  quite  good.  The  spreading 
of  the  jet  is  well  predicted  and  the  underprediction  at  the  center  line  velocity  is 
similar  to  that  observed  in  single-phase  flows  [Rodi,  1984].  Figure- 14  shows  the 
computed  v^,  at  z  s  59  mm.  The  open  circles  are  again  experimental  points,  and 

the  trend  is  similar  to  Fig.  13. 

Figure- 15  shows  the  gas  volume  fraction  as  a  function  of  the  lateral 
position.  The  open  circles  are  the  experimental  values.  The  solid  curve  is  the 
computed  Ug  at  the  axial  position  z  =  1.25  mm  (i.e.  at  the  first  gas  volume  faction 

node).  Figure  16  and  17  show  the  gas  volume  fraction  as  a  function  of  the  lateral 
position  at  distances  from  the  integration  domain  inlet  of  z  =  31  mm  and  z  =  89 
mm,  respectively.  We  can  see  that  the  agreement  is  good.  However  it  can  be 
noticed  that  the  model  tends  to  overpredict  gas  dispersion.  In  Fig.  17  the 
experimental  peaks  are  higher  than  the  predicted  ones  and  the  experimental 
center  plane  valley  is  somewhat  deeper  than  predicted. 

SUMMARY  AND  CONCLUSIONS 

A  state-of-the-art  two-fluid  model  obtained  using  ensemble  averaging  has 
been  derived  and  was  closed  using  cell  average  model.  This  approach  provides 
equations  for  multiphase  flows  that  are  mechanistically-based  (as  opposite  to 
empirical).  The  rigorous  derivation  of  the  cell  average  model  provides  exact 
constitutive  equations  for  the  invisdd  limit.  One  does  not  expect  the  values  of  the 
constants  from  cell  model  averaging  to  be  correct  for  very  viscous  flows  but  they 
provide  a  good  framework  to  start  from.  In  particular,  it  is  known  [Wang  et  al. 
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1987]  that  the  lift  coefficient,  Cl,  decreases  as  liquid  viscosity  increases.  In  this 
study,  a  lift  force  coefficient  of  Cl  =  *05  has  been  used  instead  of  the  inviscid  limit 
value  of  Cl  =  >25.  All  other  parameter  values  used  in  this  work  corresponded  to 
the  inviscid  limit  values.  The  agreement  with  the  experiments  is  remarkable. 

The  k-e  model  seems  to  be  adequate  for  this  calculation,  however  the 
observed  differences  in  phasic  dispersion  indicate  some  inadequacies  in  the 
turbulence  modeling  which  should  be  considered  in  future  studies. 
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Figure  9  Contour  Plot  of  the  Void  Fraction  (Half  of  ttie  Jet)  for  the  Vertical  Jet 

(8*0®) 


PIIOI’.NTrf 


Vector  Hot  of  the  Uquid  Velodiy  of  the 


f  the  Gas  Velocity  of  the  Indined  Jet  (0  ■  46“) 


.wWWWWWwv 


t’llOKNIC.S 


Fi2lirel2  Axial  velocity  as  a  function  of  the  lateral  coordinate  (y  -  0.1 
corresponds  to  the  jet  center  plane)  -  Inlet  conditions. 
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Void  fraction  as  a  function  of  the  lateral  coordinate 
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